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REGENERATIVE SIMULATION WITH INTERNAL CONTROLS 

by 

Donald L. Iglehart and Peter A. W. Lewis 



1 . Introduction and Summary 

Simulators are frequently faced with the task of esti- 
mating a parameter associated with the limiting distribution of 
a stable stochastic process which is being simulated. A 

methodology based on regenerative processes for obtaining point 
estimates and confidence intervals for such parameters from a 
single realization of the process has recently been developed 
in Crane and Iglehart (1974a, b), (1975a, b) and Iglehart (1975), 

(1976a, b), and (1977) . In this paper we shall introduce a new 
technique, internal control variables, which can be used in con- 
junction with the regenerative method for obtaining additional 
variance reduction for the estimates. 

Suppose regenerative process {X fc :t > 0}, which is stable 
in the sense that => X as t ■+ °° (here = denotes weak 

convergence) , is being simulated by generating a single realization 
of the process. For convenience think of this process as being 
either a positive recurrent Markov chain or the waiting time 
process for a single server queue with traffic intensity less 
than one. Then simulators are freauently interested in estimating 
r = E{f(X)}, for a given function f. The principal goal of the 
regenerative method is to produce a confidence interval estimate 
of r. The regenerative method begins by observing the pairs 



1 



is the 



of random variables { :1 _< k _< n} , where x^ 

length of the kth regenerative cycle and is the area under 

the function f (X fc ) in the kth cycle. Two basic facts are crucial 
for the regenerative method of simulating stable processes. First, 
the pairs {(Y^T^) : 1 £ k _< n) are independent and identically 
distributed (i.i.d.), and second r = E{f(X)} = E{Y^}/E{x^}. 

This last fact suggests that a natural, strongly consistent point 
estimate for r is r(n) = Y(n)/x(n), where Y(n) = n 1^=1 Y k 
and x(n) = n ^ To form a confidence interval for r 

we can use the central limit theorem (c.l.t.) 



( 1 . 1 ) 



/n (r (n) - r)/(o/E{T 1 >) = N(0,1) , 



where N(0,1) denotes a mean zero, variance one normal random variable 
2 2 

and o = E{ [Y^ - rx^] ). Of course, the constant o/E{x^} will 
have to be estimated in most simulations. 

The idea behind internal control variables is to introduce 
a third sequence of i.i.d. random variables {C^:l <_ k <_ n} which 
has the property that is defined in terms of the kth cycle 

and E{C^} is known (or can be calculated analytically) . Then 
another strongly consistent estimate for r is 



n 



-1 vn 



( 1 . 2 ) 



r CT (n) 



[Y k + 6(C k - EfC k })1 



x (n) 



The subscript CT is meant to denote "controlling the top" in the 
ratio estimator; similar estimates will be defined for "bottom control." 



2 



with a replaced 



A 

A c.l.t. analogous to (1.1) also holds for r CT 
by o', say. Since we are still free to select 8, we choose to 
do so in such a way as to minimize o ' . Having done that we find 
that 

(o ') 2 = o 2 [l - p 2 (C 1 ,Y 1 - rx 1 )] . 

To obtain significant variance reduction a control, C^ , must be 
found which is highly correlated with Y^ - ri^, a task that is 
not always easy. 

Section 2 of the paper develops the method of internal control 
variables in detail and contains specific examples associated with 
the single server queue and Markov chain models. These ideas are 
then illustrated with simulation results and numerical calcula- 
tions in Section 3. Concluding remarks are made in Section 4. 

In particular we discuss the possibility of combining internal 
controls with external controls to obtain further variance reduction. 



2 . Internal Controls 
2.1. General Ideas 

Let X = (X^_:t ^ 0} be the regenerative process being 

simulated. Recall that a process is regenerative if a renewal 

process T = (T^in 0} is defined on the same probability space 

as X and if the portions of the X process between consecutive 

regeneration points, the T 's, are i.i.d. Typically the T 's 

n n 

denote the times the X process enters a fixed state and at 
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these times the process starts over from scratch independent of 
its past history. For a formal definition of regenerative process 
and general background on the regenerative method of simulation 
analysis consult Crane and Iglehart (1975a) . Let T Q = 0 for 
simplicity and define = T R - T^, k > 1. The portion of X 

in the interval [ T k _i' T k > is referred to as the kth cycle and 
is of length x.. Suppose now that E{x.} < 00 and that the 
distribution function of x^ is aperiodic (e.g. its support 
is not contained in a set of the form { 0 , h , 2h , . . . } , where 
h > 0) Then subject to mild regularity conditions =» X 

as t Z 7 oo, where =» means P{X fc £ x} -► P(X £ x} for all 
continuity points (x) of the limit distribution. A similar result 
holds when x^ is periodic; see Crane and Iglehart (1975a). The 
random variable X is frequently thought of as the steady-state con- 
figuration of the system being simulated. Suppose f is a 
measurable function from the state space of X to the real line 
and that we are interested in estimating r = E{f(X)}. Define 
the sequence of random variables (r.v.'s) {Y^:k >_ 1} by 

T k 

(2.1) Y, = / f(X )ds , k > 1 . 

K T s — 

K-l 

If the time parameter of X is discrete (as in a discrete time 
Markov chain), the integral in (2.1) should be replaced by a sum. 

The regenerative stucture of X, the process being simulated, gives 
us the two important properties stated in Section 1: the pairs 
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{(Y k' T k ):1 - k - n * are and r = E{Y 1 >/E{t 1 } / provided 

E { I f (X) | } < °°, which we shall always assume. The c.l.t. for the 



ratio estimator r(n) = Y(n)/ T (n) indicated in (1.1) is proved 
in Crane and Iglehart (1975a) and follows from the classical 
c.l.t. for the i.i.d. mean zero, finite variance r.v.'s 

Zj, = Y^ - ri^, k >_ 1 . We always assume 0 < = E{Z^} < <». 

2 ^ 

The variance of Z ± , o , is also related to the variance of r(n) 

through the asymptotic relation (as n -> °°) 



(2.2) o 2 {r (n) } - n 1 (o 2 {Z 1 }/E 2 {t 1 } ) . 

For a derivation of this result see Cramer (1964, p. 354, 
eq. (27.7.3)). A number of point estimates and confidence intervals 
have been proposed for ratios such as r (see Iglehart, 1975) , but 
the simplest conceptually and computationally are the so-called 

/*s 

classical ones. The point estimate is r(n) and the confidence 
interval is 



(2.3) 



I(n) = [? (n) - z i_ y / 2 S//j ' r < n ) + 2 1 - y / 2 S//t -* ' 



where z, = $ ± (1 - v/2) , the inverse of the standard normal 

l-y/2 

distribution function, and s is the classical point estimate 
of o which is constructed as follows. Let s^ [resp. s 09 ] 



22 



A 

be the sample variance of the ’ s [resp. t^'s] and s^ 2 

the sample covariance of the Y^ 1 s and t^'s. Then s is defined 

by 

s = [s 11 - 2 (Y/t)s 12 + (Y/t) 2 s 22 ] 1/2 . 
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Introduction of the internal control variables mentioned 

in Section 1 is done with the hope of reducing the variance of 

the point estimator of r. This in turn, will either reduce the 

number of cycles that need to be simulated for fixed precision 

or reduce the length of the confidence interval, (2.3), for r. 

The sequence of internal control variables {C^:k 1} is defined 

on the same probability space supporting the process X in such 

a way that only depends on X (or underlying r.v.'s 

defining X) in the kth cycle. This construction of the C, ' s 

will insure that they are i.i.d. because of the basic i.i.d. 

structure of the cycles of a regenerative process. The control 

variables C, allowed, however, are further restricted to those 
k 

for which the mean, EtC^}, is either known or can be calculated 

analytically. Thus one of the prices we must be prepared to 

pay to obtain variance reduction for our simulation is the 

analytical work of computing E{C^} . Having defined the ' s 

we proceed to form new ratio estimators for r. Starting with 

the ratio estimator r(n) = Y(n)/T(n) we can either 

control the Y, 's, the t 's, or both if we introduce a second 
K k 

sequence of control variables. The estimator r CT (n), proposed 
in (1.2), involves controlling the Y^ ' s only. Here the sub- 
script CT is meant to suggest "controlling the top." For 
convenience let (for k >_ 1) 

Y k - \ + 8(c k - E(c k }) 

and 
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Observe that the 1 s are i.i.d. with E{Y|) = E{Y^} and con- 
sequently the Z£ ' s are also i.i.d. with E{Z£) = 0. Hence the 
standard c.l.t. yields (as n -*■ °°) 

l Z '/o{Z' }n 1/2 = n 1/2 [ (Y/t ) - r] / (a{ Z ' }/x) =» N(0,1) 
k=l K 1 1 

which, upon replacing x by E{x^} in the denominator (this 
can be justified by a continuous mapping argument) , oecomes 

n [r (n) - r] 

(a{ Z| j/Elx'jTT “ N(0 ' 1) * 

The variance of Z^, which obviously depends on B, is easily 
seen to be 

cr 2 { Z|) = cr 2 {Z^} + 2B cov{Z^,C^} + B^a 2 {C^} . 

2 

Now select B so as to minimize a {Z'K This yields 

cov{Z , C. ) 

(2.4) B* = - --5 — 

a {C^} 

and 

(2.5) o 2 { Zp = a 2 {Z 1 > [1 - p 2 {Z 1 ,C 1 )] , 

where p{Z^,C^} is the correlation coefficient between Z^ 

and C 1 . If we use the c.l.t. for the estimator r CT (n) as a basis for 



7 



constructing confidence intervals for r, then we want a control 

2 

variable which will minimize a {Z|}. From (2.5) we see 

2 

that we need a large value of p {Z^C^}. Since the length 
of such a confidence interval is proportional to a{ Z| } , 
to reduce the number of cycles required by a factor of four we 
need |p{Z^,C^}| > (0.75)^^^ = 0.866. To obtain a control, C^, 
which is highly correlated with either or is relatively 

easy to do. However, because and are themselves highly 

correlated, it is much more difficult to find a control 
which is highly correlated with Z^ . In general, we shall try 
to find a control which mimics Z^ but for which we can 

still calculate its mean, Elc^) . As usual we try to do as much 
analytically as possible, before having to resort to simulation. 

We illustrate some possible candidates for controls in the context 
of GI/G/1 queues and discrete time Markov chains. 



2.2. Internal controls for GI/G/1 Queues 

In the GI/G/1 queue we assume the zeroth customer arrives 

at time t^ = 0, finds a free server, and experiences a service 

time Vq . The nth customer arrives at time t and experiences 

a service time v . Let the interarrival times t - t . = u , 

n n n-1 n 

n > 1. Assume the two sequences (v :n > 0} and {u :n > 1} 

— n — n — 

each consists of independent, identically distributed (i.i.d.) 
random variables (r.v.'s) and are themselves independent. Let 
E{v n ) = p 1, E{u } = A "*■, and p = A/y, where 0 < A , y < 00 . Thus 
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y has the interpretation of the mean service rate and X has 

the interpretation of the mean interarrival rate. The parameter 

p is called the traffic intensity and is the natural measure 

of congestion for this system. We shall assume that p < 1, a 

necessary and sufficient condition for the system to be stable. 

While many characteristics of interest can be estimated 

using the regenerative method, we shall restrict our attention 

th 

to the waiting time of the n — customer, W n (time from arrival 

to commencement of service) . For further discussion of the 

simulation of the GI/G/1 queue see Crane and Iglehart (1974b) . 

To obtain a representation for the process {W^in >.0} let 

X = v . - u and set = 0, S =X.+***+X,n>l. 

n n-1 n 0 n 1 n — 

The following well-known recursive relationship exists for 

the W ' s : 
n 



W n = 0, 



w i — [w 

n+1 n 



X n+1 ] ' 



n > 0 



By induction, we also have 



W = max(S - S. :k = 0, 1, ... , n}, n > 0. 

n n k — 

Using the strong Markov property of the process ^ s n :n 2. ^ } 

it can be shown that there exists a sequence of integer-valued 

r.v.’s (T^:k 2. 0} such that T^ = 0 , T^ < T^ + ^, and W T = 0 

k 

with probability one. In other words, the customers numbered T^ 
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are those lucky fellows who arrive to find a free server and 

experience no waiting in the queue. The fact that there exists 

an infinite number of such customers in the GI/G/1 queue is a 

direct consequence of the assumption that p < 1 and the 

strong law of large numbers. Thus {W^in ^ 0} is a regenerative 

process. If we let T k = T k ~ T k-1' ^ 1. b • then represents 

th 

the number of customers served in the k— busy period (b.p.) 
and they are numbered T k-1 + • Next define 

the sequence {Y^:k 'L 1) by 



Y 



k 



Tk ' 1 

l 

3 =T k-l 



W, 



t 



the sum of the waiting times in the k— b.p. Since the queue 

is stable for p < 1, we know that W =» W as n -► 00 and wish 

to estimate E{W} . 

We are interested in constructing controls, , which are 

highly correlated with z i =Y i~ rT i' but which are not so complex 
that we can not calculate their means. The controls we consider 

are of the form c i = D q - where the r.v. is an 

attempt to mimic . To compute the mean of we need, of 

course, to be able to compute E{t^}. For M/G/l queues we 
know that E{t^} = l/(l-p). For GI/M/1 queues E{t^} = l/(l-6), 

where 6 is the root inside the unit circle of z-E{exp [-y (1-z) u^ ]} = 0 
where u^ is an interarrival time . If the queue in question is 
neither an M/G/l or GI/M/1 queue, then the term i^/y in 
will have to be approximated by another r.v. whose mean can 
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be computed. The factor 1/p multiplying x^ helps to make 
the variance reduction obtained independent of the scale parameter 
p. Also it can be argued that the term x^/p is then in the 
same units as D^, namely, the unit of time. Several alternatives 
for are listed below, indexed by a superscript. They are 



d' 1 ’ - x; - 1 



w o = 0 ' 



w o + w i ' 



T i = 1 



T 1 > 2 ; 



D 



( 2 ) 



X 1 + 4 ’ 



T 1 = 1 



T 1 1 2 ' 






W 0 = 0 ' 

w 0 + w x , 



w 0 + W 1 + X 2 ' 



T 1 = 1 
T 1 = 2 



r 1 > 3 ; 



and 



°1 3) = (X 1 + X 2 )+ = W 2 ; 



D 



( 4 ) 



0 , 



X+ + (X+ + x 2 ) + 



T 1 = 1 



x 0 >_ 2 






W 0 = 0 ' 

w 0 + w x , 

w 0 + w x + W 2 , 



T 1 = 1 
T 1 = 2 
T 1 > 3 
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In general, it is more difficult to calculate E{D^^} as i 
increases. On the other hand, as i increases comes 

closer to and presumably results in more variance reduction. 

Let = D^ 1 ^ - T^/y, i = 1,2,3. 

In our simulation runs for an M/M/1 queue the controls 
, C< 3) , and were generally found to do much better than 

. However, the more complicated controls, and , 

( 2 ) 

gave little improvement over . Thus with both variance 

( 2 ) 

reduction and ease of computation in mind, is our first 

( 2 ) 

choice. To compute for the M/M/1 queue first note that 

oo 

P{t = 1} = p{v < u, } = J e ^ ye dy = (1 + p) 1 , 
i u i Q 



and 



Hence 



<t£x> f d Y 



E{X 1 |X 1 > 0} W U ' 



E{X+} = E{X 1 IX 1 > 0}- p{X 1 > 0) = ;nT £ ?T 



e{d 




O + t-2- e(x+ + X+|x* > 0} 

±+p 1 Z 1 



rl + p i 

1+p L p + p(l+p) J ' 



( 2 ) 

since X^ and X 2 are independent. The expectation of C-^ 
is then given by 
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2 



E{cj 2) }= u 1 



[ 



1+P 



^i+p ; 



-&] ■ 



Further discussion on the computation of the mean of controls such 
( 2 ) 

as for the GI/G/1 queue is contained in the Appendix. Com- 

putational results from our simulation runs are contained in 
Section 3. 



2.3. Internal controls for Markov chains 

The second example we consider is a discrete time Markov 
chain. Assume now that we are simulating an irreducible, aperiodic, 
positive recurrent Markov chain, X = {X n :n 0}, with the goal 
of estimating r = E{f(X)}. The controls we consider here are 
of the form 

n o A ( x i - D 

(2.6) C = l [f (X £ ) - r ] , 

1 £=0 0 

where r^ is a guess of r, n^ a fixed integer, and 

n 0 A denotes the minimum of n^ and ( T ^-l) . Again our 

motivation for is to mimic 

T i _1 

z = l [f <X ) - r) = Y 1 - rt 1 

1=0 

Of course, we still must be able to compute E{C^} in order 
to implement the method of internal control variables . Let ^ 
be the state space of X, which we assume is finite, and P the 
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transition probability matrix. Furthermore, let g be the 

column vector with components f(i) - r„ and P be the matrix 

0 3 

with elements 



( Pk£' 

jPk£ = j 

' 0 , 



* * j 

£ = j . 



Then if we let = j and form regenerative cycles based on 

the times of return to state j, it can easily be shown, using the 
method developed in Hordijk, Iglehart, and Schassberger (1976), 
that 

n o 

(2.7) E(C } = ( l P n g) , 

n=0 3 3 



where j is the regenerative state and the subscript j means the 
jth component of the indicated vector. For two simple Markov chains 
the theoretical amount of variance reduction obtained using the 
control has been calculated for different regenerative states j 

and values n^ and is contained in the next section. From (2.6) 
it is clear that the larger the integer n^ selected the closer 
C-^ comes to duplicating . However, the larger n^ the more 
difficult is the computation of E(C^) contained in (2.7). 
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3. Numerical Results 



3.1. The M/M/1 queue 

Even for the simple M/M/1 queue it is in general impossible 
to verify analytically what variance reduction will be obtained 
via the several internal controls listed in the previous section, 
or to get an idea of the magnitude of the effect. For something 
as simple as it is difficult to compute analytically the 

correlation between and for the M/M/1 queue, and this 

is what is required in equation (2.5) to find the variance 
reduction. 

Thus, we resorted to simulations to verify the amount of 
variance reduction obtained and the relative effectiveness of 
the various controls. In the final simulations all runs were 
performed on an IBM System 360/67 computer using the LLRANDOM 
package (Learmonth and Lewis (1973)) which generates random numbers 
according to the scheme given by Lewis, Goodman, and Miller (1969) 
and exponentially distributed random numbers using the Marsaglia 
"rectangle-wedge- tail" method. Tests of the random number 
generator are given in Learmonth and Lewis (1974). 

Of the extensive simulation runs performed, we give here 
only a summary of the conclusions and two detailed tabulations 
in the case of the most suitable control. 

(1) The controls and ^ do much better generally 

than D {"^ ' with little improvement over obtained by 

use of and D i^ * We sa y generally because results 

vary unpredic tably with \ and y and their ratio o. 
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(2) Subtracting the number of customers served in a busy period 

generally improves the variance reduction. By making it 

( 2 ) 

dimensionally stable as in we obtain a "variance 

reduction" measured in terms of ratios of standard deviations, 

of approximately 70%, uniformly over A and y. This is 

roughly equivalent to halving the number of cycles (b.p.'s) 

2 

that one must simulate; (0.7) ~ .5. Much better reductions 

can be obtained for smaller p (i.e. p = 0.25) by specially 

( 2 ) 

designed controls; the point is that works even out 

at p = 0.99, where variance reduction is extremely important. 

Table 1 shows results obtained by simulating an M/M/1 

( 2 ) 

queue out to n = 2000 cycles with control C 1 and replicating 
the simulation either m = 250 or m = 100 times to estimate the 
variance of the estimators r(n), r^,p(n), and r^g(n), where we 
drop the n for convenience. Here, we have specifically that 



n 



r CB (n) “ 



n k=l K 



1 n 

- y 

n t* _ 



n ■ [T k + S(C 1 2> 
11 k = l K 1 



E{c| 2) }) ] 



The estimated precision (standard deviations) of the estimates 
of E(W) are given in brackets under the estimates. 

The results in Table 1 are for p = 0.5 (and m = 250) 
and two values of y and for p = 0.99 (and m = 100) 
for two values of y commensurate with those given for 
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p = 0.5. The second, third and fourth columns in the Table give 
estimates of the correlations between the control and Y^-t^ etc. from 
which the theoretical variance reduction can be computed. They 
are very close to the values given in the next to last column, 
the observed variance reduction, from which we deduce that 
estimating 3^ and 8g affects the variance reduction only 
slightly. Overall there is negligible effect of different values 
of y on the variance reduction obtained for the cases p =0.5 
and p = 0.99. For the results p = 0.99 given in Table 1, the 
variance reduction is 75%, which is about the same as for 
p = 0.5. For the case where the control is on the bottom, i.e. 
for x, the variance reduction is not quite as good for control 
of Y. Note too that the estimated values of E{w} appear in 
some cases to be at least three or four standard deviations from 

A A A 

the true mean. This is because the estimates r, r and r 

Cl Ld 

can be seen from the 100 replications to be non-normal. In other 
words, for high p(0.99), the simulation needs to be taken out 
further than 2000 cycles. In summary, the variance reduction 
obtained with the internal control is practically independent 
of the scale factor y and the traffic intensity p. 

3.2. Two simple Markov chains 

We turn now to two simple Markov chains to qive another 
illustration of the method of internal controls. The first 
chain is simply a recurrent random walk with state space 
E = {0,1, 2, 3, 4} and transition matrix 
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1/2 


1/2 


0 


0 


0 


1/4 


1/2 


1/4 


0 


0 


O' 


1/4 


1/2 


1/4 


0 


0 


0 


1/4 


1/2 


1/4 
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Suppose we are simulating to estimate r = E{ X} , which can be 
computed to have the value 2. Table 2 contains the theoretical 
values of a {Z|}/a {Z^} and p(C^,Z^), where is given by 

(2.6) for several choices of n^ and r^ . The entries in these 
tables were calculated using the methods in Hordijk, Iglehart, 
and Schassberger (1976) . Note that the variance reduction is 
largest when E{t^} is smallest and n Q is largest. This also 
yields the highest value of p(C^,Z^) in keeping with our idea 
of having mimic Z^ . Not too much is lost in variance 

reduction as the guess, r^, departs from the true value r = 2. 
A variance reduction of more than 50% is attainable, but depends 
on the state which is used as the regeneration point. 

The second Markov chain represents an (s,S) inventory 
model; see Crane and Iglehart (1974) for further details. The 
state space $ = {6,7, ...,10} and the P matrix is given by 
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TABLE 2 



2 2 

Variance ratio o Z| /a {Z^} (top row) 



and 


Correlation 


P ( C j ) 


(bottom row) 
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TABLE 3 



2 2 

Variance ratio a {Z /a { Z^ (top row) 
and Correlation p(C^,Z^) (bottom row) 
for (s,S) Inventory Model 
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Again we wish to estimate E{X} = 8.297. Table 3 contains the 
results analogous to those for the random walk (Table 2) . The 
same general observations about E{t^} and n^ made for 
Table 2 seem to hold here as well. Note however that the 
results are sensitive to the choice of r^ . Again a 50% 
reduction in variance can be attained. 



4 . Conclusions and Extensions 

We have been able to obtain a 50% variance reduction using 
internal control variables, for the regenerative estimate of the 
limiting value of the mean waiting time in an M/M/1 queue. This 
reduction is obtained uniformly over all parameter values. It is 
fairly certain that the technique will work well with any GI/G/1 
queue . 

Internal control variables can be easily used with 
discrete time Markov chains. The examples used in this paper 
showed that a variance reduction of 50% is attainable. This figure 
is likely to vary widely with the particular Markov chain. 
Continuous time Markov chains and semi-Markov processes can be 
handled in the same way using the discrete time method of 
Ilordijk, Iglehart and Schassberger (1976) . 

Another method of internal stratified sampling was also 
investigated. This method produced little overall variance 
reduction despite considerable effort. 
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Finally we note that it is possible to apply the idea of 
internal controls to the classical sample average estimate over 
a realization of fixed length m. Thus in estimating E(W) 
in the GI/G/1 queue we have 

, m-1 

W (m) = — l W. , (4.1) 

m j=0 11 



which may be written as 



N(m) 

W(m) = m { 2 Y k + Y - 

m k=l K 



n(m)+l } ' 



(4.2) 



where N(m) is the number of completed busy periods in the queue 
in [0,m], as before is the sum of the waiting times in the 

kth cycle, and Y' . is the sum of the waiting times in the 
last, incomplete cycle. 

A central limit theorem for W (m) holds for which the variance 

2 

term is proportional to E{Z^}, and which is now estimated 
from the random number N (m) of cycles. Control is applied 
to the Yj, ' s in (4.2) just as it is applied in the ratio 



estimator r^ T (n) . Call this estimator W^ (m) . For the M/M/1 
queue the variance reduction observed in the simulations was 
the same for all controls C^^ with the ratio estimator and 
with W c (m) . 



The main reason for considering (m) is that, while 
it loses the advantage of being an estimator using a fixed number 
of i.i.d. random variables, one can apply the classical external 
controls to (m) . Thus one could use the difference of the sum 
of the m arrival times and the m service times to control 
W (m) . We have not yet tried this idea. 



23 



APPENDIX 



To implement the internal control techniques for a 
GI/G/1 queue, certain theoretical parameters are required. In 
this appendix we shall indicate the values of these parameters 
in so far as they can be calculated. These values are either 
well-known or easily calculated. For a reference to the known 
formulas see Cohen (1969) . 

We begin with E{x^}, the expected number of customers 

served in a busy period. For the general GI/G/1 queue recall 

that we let X = v - u and S = X, + • • • + X , for n > 1, 
n n-1 n n 1 n — 

with S rt = 0 . Then x, = inf{n > 0:S < 0}. The general 

0 1 n — 

expression for E{x^} is given by 

oo 

E{x 1 > = exp{ l n -1 P[S n > 0]} , 

an impossible expression to evaluate in general. Another useful 
expression for E{x^} is 

(A . 1 ) E { x x } = 1/P{W = 0} , 

where W is the stationary waiting time. In the special case 
of the M/G/l queue, however, we have 

E{ Tl } = (1 - p)" 1 . 
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Now for the queue GI/M/1 we can use (A.l) and the stationary 
distribution of the embedded Markov chain to conclude that 

Et-^} = (1 - 6 ) -1 , 

where 6 is the root inside the unit circle of 

z - U{y (1-z) } = 0 

-su 

with U(s) = E{e }, Re s > 0, and where Vq is an exponential ( y ) 
r.v. It is easy to check that 6 = p for M/M/1 queues. Daley 
(1975) has recently proposed the approximation to 6 given by 

6 = a 1 ( 1 — p ) 2 + 2 (l-b _1 ) p + (2b _1 -l) p 2 , 

2 

where a^ = P{u^=0}, E{u^} = 1, and b = E{u^} . This approximation 
gives good results in a number of examples calculated by Daley 
(1975) and may be useful for the purposes of this paper. 

Next we turn to the computation of P(x^=l} and P{t^=2) 

For the GI/G/1 case we have 

P{t 1 =1} = P^ £ 0} 

and 

P{t 1 = 2} = P{S 1 > 0, S 2 <_ 0}, 
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both of which can be worked out with a little effort. For the 
M/M/1 queue 

P{t 1 =1) = (1+p) ^ 

and 

P{t 1 =2} = p(l+p) -3 . 

For the M/G/l queue 

P{x 1 =l) = V ( A ) , 

_Xv o 

where V(A) = E(e ), while for the GI/M/I queue 

P{t 1 =1} = 1 = U (y ) , 

where U(s) is given above. For the M/E^/l and E^/M/l queues 

the value P{x^=2} can be calculated with some effort. As these 

expressions are cumbersome they shall be omitted. 

Finally we give various partial expectations which are 

( 4 ) 

needed for computing the means of internal control . Namely, 

E{S+ + S+, _> 2} = E{S lf S 1 > 0} + EfS*, S 1 > 0} 

and 

E{x l' t i 1 2} = - P{x 1 =l}. 

Here the symbol E{X,A} = E(Xl a ), where X is a r.v., A an 
event, and 1^ the indicator function of A. In the special 
case of the M/M/1 queue. 
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E{S X , S > 0} = p/{p(l+p)} , 



E(S 2' S 1 > 0) = [ 2 (if?) 



2 2 
, P 

1 

(1 + p 



] " • 



and 



E ^ T 1' T l — = 2p/{ (1-p) (1+p)} . 
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